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Abstract 

The existence of analytical models that foresee the pressure behavior during injectivity tests 
is required to estimate reservoir parameters and properly interpret pressure transient data. 
However, the currently available formulations can only describe injectivity tests in single¬ 
layer reservoirs with horizontal wells and multilayer systems with vertical wells. Therefore, 
the main purpose of this work is to provide a model that depicts the pressure behavior of 
multilayer reservoirs with multilateral horizontal wells under two-phase flow. A new real 
space solution for single-phase flow in those systems was also developed. The proposed 
formulation was also used to estimate the reservoir equivalent permeability. 


1 Introduction 

Adequate reservoir resources management relies on knowledge regarding its parameters, 
such as permeability and extension. One conventional way to estimate those parameters 
consists of well testing, during which the reservoir underlying fluid is produced and pressure 
data is measured. At a first stage, the wellbore is open, the reservoir underlying fluid 
is produced and, thus, pressure decreases. Afterwards, the wellbore is shut and pressure 
rises (Gao, 1987; Goode & Thambynayagam, 1987). Those stages are called drawdown and 
buildup, respectively (Ehlig-Economides & Joseph, 1987; Odeh & Babu, 1990). 

Injectivity tests are an alternative to conventional well testing. This procedure is based 
on pressure transient analysis due to the injection of an external fluid into the reservoir 
(Peres & Reynolds, 2003). Injectivity tests present some advantages compared to well 
testing, such as the dismiss of the produced fluid disposal, reduced environmental fingerprint 
and increased operational safety (Bonafe, 2019; Salina Borello et ah, 2019). Analogously to 
well testing, two distinct stages are also observed at an injectivity test: injection and falloff. 
During the former, a fluid is injected into the reservoir and pressure rises. The latter is 
identified by a pressure drop, after the well is shut (Peres et ah, 2004; Bela et ah, 2019b). 

Injectivity tests may be performed either in vertical or horizontal wells. The latter 
present some advantages compared to the former, such as: reduced fluid coning and better 
performance in reservoirs with high vertical permeability (Ozkan et al., 1989). Moreover, 
horizontal wells can also enables drilling operations underneath a permanent structure, such 
as an airport, and a single horizontal wellbore can replace several vertical wells (Zhan & 
Zlotnik, 2002). 

Reservoir parameters are estimated from analytical models that depict the pressure 
behavior during the test. In multilayer reservoirs with vertical wells, the solution for both 
single (Lefkovits et al., 1961; Gao, 1987; Ehlig-Economides & Joseph, 1987) and two-phase 
flow (Peres et al., 2004; Barreto Jr. et al., 2011; Bela et al., 2019b; Bonafe, 2019) is already 
known. For horizontal wells, on the other hand, currently existing formulations are restricted 
to single-layer reservoirs (Goode & Thambynayagam, 1987; Daviau et al., 1988; Ozkan et 
al., 1989; Peres & Reynolds, 2003; Peres et al., 2004). The available solutions for multilayer 
systems with horizontal wells are limited to single-phase flow and were developed in Laplace 
domain (Kuchuk & Habashy, 1996; Larsen, 2000). 

Hence, this work aims to propose an analytical model for the injection period of injec¬ 
tivity tests in multilayer reservoirs with horizontal wellbores. Since there is currently no 
real space formulation for fluid production through multilateral horizontal wells in stratified 
systems, first a single-phase solution was reached using source functions and Newmann’s 
product. Then, the proposed analytical model was extended for two-phase flow as well, by 
adapting the single-layer solution achieved by Peres and Reynolds (2003). An adjustment 
is also suggest to apply the average pressure technique in multilayer systems. The proposed 
model was validated via comparison with a commercial finite difference-based flow simu¬ 
lator. Besides, the reservoir equivalent permeability was estimated through the pressure 
derivative analysis. 
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A brief overview about previous works regarding analytical solutions for horizontal wells 
is conducted in section 2. Section 3 presents the suggested single-phase flow formulation, 
while the two-phase model is depicted in section 4. The proposed method was applied on 
a set of cases, whose results are presented in section 5. Last, the main conclusions of this 
work are reported in section 6. 


2 Previous Achievements 

The first analytical models for single-phase flow through horizontal wellbores were pre¬ 
sented in the late 80’s. Goode and Thambynayagam (1987) achieved a solution for pressure 
behavior through horizontal wells by applying successive Laplace and Fourier transforms. 
They also proposed approximate expressions for each observed flow regime. Their formula¬ 
tion applies the superposition principle to compute buildup pressure. 

Daviau et al. (1988) proposed a solution for single-layer reservoirs based on Green 
functions and the sources/sinks method. The formulation assumes uniform influx wellbore 
condition, that is, fluid inlet along the wellbore is constant. They also made brief comments 
regarding the different flow regimes observed and the influence of wellbore hydraulics and 
positioning. 

Following their work, Ozkan et al. (1989) provided further insights over the existence 
of two distinct radial flow regimes that may be detected at well testing with horizontal 
wells. Anisotropy was included and they suggested a technique to estimate the reservoir 
permeability and skin factor from the pressure derivative with respect to the natural loga¬ 
rithm of time. Besides, it was stated that the pressure response at a wellbore with infinite 
conductivity may be estimated by the equivalent pressure point of an analytical model that 
assumes uniform flux. 

An analytical model for the infinite conductivity was developed by Rosa and Carvalho 
(1989). Their hypothesis considers that pressure in the wellbore interior is constant. To 
solve the problem, they proposed that the well may be divided into n segments, each of 
them presenting uniform influx. The pressure equality in all segments, combined with the 
constant total producing flow-rate, yield a linear system that allows the determination of 
flow-rate in each segment. 

Approximations for each flow period and the required conditions for them to develop 
were summarized by Odeh and Babu (1990), as well as estimates for their starting and ending 
times. Once again, buildup pressure is computed through the superposition principle. 

The solution for horizontal wells in Laplace domain was reached by Kuchuk et al. (1991). 
Their model is flexible when it comes to the reservoir outer boundaries conditions, since it 
foresees cases with constant pressure maintenance systems, such as gas caps, acting at one 
of the reservoir vertical boundaries. 

Jelmert and Thompson (1991) proposed that the linear system associated with the 
infinite conductivity condition is easier to compute in Laplace domain than in the real 
field. They studied the influence of assuming a piecewise linear flow-rate distribution, when 
compared to a stepwise constant flow-rate distribution, as suggested by Rosa and Carvalho 
(1989). Their results supported that it is possible to apply the equivalent pressure point of 
an uniform influx model to compute the pressure response of an infinite conductivity well. 

A generalized solution for slanted wells was reached by Zhan and Zlotnik (2002). Their 
solution consists of dividing the wellbore into several sink points. Then, the solution for those 
points is determined in Laplace domain. Horizontal and vertical wells may be expressed as 
particular cases of their formulation. 

Wang and Zhan (2017) have evaluated the influence of wellbore hydraulics into the 
pressure response. They compared the equipotential curves in a given reservoir assuming 
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three distinct wellbore conditions: uniform influx, infinite conductivity and mixed boundary 
condition, which describes the wellbore hydraulics in a more realistic way. They showed that 
although pressure change behaves as a linear function of flow-rate for both uniform influx 
and infinite conductivity, the same does not hold for a mixed boundary condition. 

When it comes to multilayer reservoirs, a formulation in Laplace domain was developed 
by Kuchuk and Habashy (1996). Their solution assumes that the wellbore perforates only 
one layer and, as the test goes on, pressure behavior gets progressively influenced by the 
other layers, due to the formation crossflow. 

An analytical model for multibranched horizontal wells was proposed by Larsen (2000). 
This solution, which was also reached using Laplace transform, approximates each wellbore 
ramification by a uniform influx fracture and is applicable in multilayer reservoirs, with or 
without formation crossflow. 

As far as the authors can state, the formulations achieved by Kuchuk and Habashy 
(1996) and Larsen (2000) are the only existing solutions for single-phase flow in multilay¬ 
ered reservoirs with horizontal wells. However, both models were developed in Laplace do¬ 
main. Therefore, a numerical inversion algorithm is required to apply their solutions. Thus, 
the first goal of this work consists of extending the real space formulation for single-layer 
reservoirs (Daviau et al., 1988; Ozkan et al., 1989) to multilayered systems. 

Only a few papers regarding two-phase flow in horizontal wells have been published. 
Peres and Reynolds (2003) developed a formulation that describes the pressure response 
during the injection period for single-layer reservoirs. They proposed that the mobility 
differences between the fluids may be condensed into a term that adjusts the expected 
single-phase flow pressure change. This additional term is computed after the flood front 
in each direction has been determined. Their results also showed that pressure response 
from fluid injection through a horizontal wellbore might be estimated by the single-phase 
solution, unlike what happens in vertical wells. 

The falloff solution in single-layer reservoirs was reached by Peres et al. (2004), applying 
the superposition principle to depict flow-rate history after the well is shut. Both injection 
and falloff solutions assume uniform influx along the well. 

To the best of our knowledge, there is currently no analytical solution for the pressure 
behavior in multilayer reservoirs under two-phase flow. Thereby, the main purpose of this 
work is to achieve a model for the injection period during injectivity tests in stratified reser¬ 
voirs with multilateral horizontal wells. The proposed formulation is based on the already 
known single-layer solutions (Daviau et al., 1988; Ozkan et al., 1989; Peres & Reynolds, 
2003; Peres et al., 2004). 


3 Single-Phase Flow Through Multilateral Horizontal Wells 

The suggested formulation considers a reservoir with an arbitrary number of layers, 
which may present distinct properties. It was assumed no formation crossflow. Therefore, 
pressure behavior will be influenced by all layers only if the wellbore presents one ramification 
per layer. Figure 1 illustrates the considered reservoir model and wellbore configuration. 
All computations assume that a consistent set of units is used and the following simplifying 
hypothesis: 


• Reservoir initially in equilibrium; 

• Homogeneous, laterally infinite reservoir, with impermeable vertical boundaries; 

• Multilayer reservoir without crossflow, that is, the wellbore is the only flow path 
between layers; 

• Apart from the hydrostatic column, pressure is initially the same in all layers; 
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• Wellbore is parallel to the horizontal plane, located at an arbitrary distance from the 
formation bottom and produces at a constant flow-rate q ; 

• Infinite conductivity: pressure is constant along the well; 

• Fluids are slightly compressible and present constant viscosity /z; 

• Isothermal flow; 

• Gravitational and wellbore storage effects are neglected; 

This section briefly presents the solution for single-phase flow in single-layer reservoirs 
with horizontal wells. Then, it will be extended to multilayer systems. 

3.1 Single-phase flow in single-layer reservoirs 

The solution for single-phase flow in single layer reservoirs is well known (Daviau et ah, 
1988; Ozkan et ah, 1989). Considering that the wellbore heel is positioned at coordinates 
(xi,yi, z\) and applying Newmann’s product, the pressure change A P at a given time t may 
be decomposed as the product of three terms, such that each term is associated to flow in 
a given direction (Rosa & Carvalho, 1989): 



o 


In equation (1), q stands for the flow-rate and the wellbore length is represented by 
L 1 while <f> and c* denote the porosity and total compressibility, respectively. Equation 
(1) represents the time integral of all instant pressure variations at a given point with 
coordinates (x,y,z) due to the source point productions along the well. Instant pressure 
change in x-direction are given by (Daviau et ah, 1988; Ozkan et ah, 1989): 



( 2 ) 


where y x = k x /ifi/ict stands for the hydraulic diffusivity coefficient in x-direction. Anal¬ 
ogous hydraulic diffusivity coefficients are also defined in y and z-directions. 

In y-direction, the contribution of all point sources along the wellbore must be accounted 
for. Thus, considering that the wellbore section open to flow goes from y = y\ to y = yi +L: 



z 



Figure 1 . Multilayer reservoir model with a multilateral horizontal wellbore 
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Vertically, the reservoir is finite. Therefore, the images method must be applied in 
z —direction (Rosa & Carvalho, 1989). Thus, denoting the reservoir thickness as h and the 
smallest distance between the wellbore and a vertical boundary as dz: 


A P z {z,t) 


1 

2y/n'q z t 


OO 

^2 exp 

.71 —— OO 


(2 nh) 2 \ 

4 ?ut J 


+ exp 


/ (2nh — 2dz) 2 \ 

V 4r? 2 f ) 


( 4 ) 


Equations (1) to (4) are achieved assuming uniform influx, that is, flow-rate per unit 
length is constant along the well. However, flow in horizontal wells is more accurately 
depicted by the infinite conductivity hypothesis, which states that pressure is constant 
along the wellbore (instead of fluid in flow). Pressure behavior in an infinite conductivity 
well may be estimated from a uniform influx model, using the average pressure technique 
or the equivalent pressure point method (Goode & Thambynayagam, 1987; Daviau et al., 
1988; Ozkan et al., 1989; Kuchuk et al., 1991). 


3.2 Single-phase flow in multilayer reservoirs 

The formulation for multilayer systems starts by applying equation (1) into a given 
layer j: 


A Pj(x,y,z,t) 



t 

J A P Xj (x,T)AP yj (y 1 T)AP Zj (z,T)dT 
0 


( 5 ) 


where the directional pressure changes are computed analogously to equations (2) to (4). 

Equation (5) assumes that layer flow-rate is constant in time. Otherwise, layer flow-rate 
qj would be a function of time and could not be written out of the integral sign. Rigor¬ 
ously speaking, if layer properties are remarkably different, pressure diffusion propagates 
at differently in each layer and a flow-rate transient may develop. However, based on the 
formulation for multilayer reservoirs with vertical wellbores (Gao, 1987; Ehlig-Economides 
& Joseph, 1987), this rate transient period is not expected to last long. The influence of 
this hypothesis will be further analyzed in section 5.1. Therefore, condensing the integral 
term in equation (5) into a new auxiliary parameter <jj(xpy, z,t), flow-rate in layer j may 
be written as: 


= Ap ( x y z t y 

aj(x,y,z,t) 


Uj (x, y, z,t) = / A P x AP y . AP z dr (6) 


By model hypothesis, pressure change is the same in all layers (except for the hydrostatic 
column). Thus, summing all layer flow-rates: 


LjcpjCtj 
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( 7 ) 


Rearranging equation (7) in a more convenient way: 
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AP(x,y,z,t) = q 



LjipjCtj 

aj(x,y,z,t) 


-1 


( 8 ) 


For single-layer reservoirs, equation (8) reduces to equation (1), assuring that the pro¬ 
posed formulation is applicable to reservoirs with any number of layers. Finally, to account 
for the skin effect, an additional pressure change must be added. In horizontal wells, the 
damage region affects flow in zx plane. For this reason, the mechanical skin is computed in 
terms of the wellbore total length (Daviau et al., 1988; Kuchuk et ah, 1991): 


AP(x,y,z,t) = q 



Lj <t>j ctj 

aj(x,y,z,t) 


-1 


qy-S . 

(kzx')eql't 


E QjSj 

S = — - (9) 
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In equation (9), y represents the reservoir underlying fluid viscosity and L t denotes 
the sum of all wellbore branches lengths, while S stands for the reservoir overall skin factor 
(Gao, 1987). The overall skin factor accounts for all layer skin factors Sj (as defined by 
Hawkins (1956)). The equivalent permeability in zx plane is denoted by ( k zx ) eq and is 
defined as a length-averaged permeability: 


n 

E Lj \Jk Xj k Zj 

(k zx ) eq = — - j - ( 10 ) 

3.3 Comments about some simplifying assumptions 

Equation (9) may be used to compute the pressure change at a given point of a reser¬ 
voir subject to the assumptions mentioned above. No restriction is made regarding the well 
positioning. That is, the developed formulation accounts for wellbore off-centering. Proper¬ 
ties such as permeability, porosity, skin factor and well length can also be different for each 
layer. The strongest hypotheses made to achieve equation (9) consist of assuming uniform 
influx and constant layer flow-rates. 

Nevertheless, the proposed formulation can easily represent an infinite conductivity 
wellbore using either the average pressure technique (Goode & Tliambynayagam, 1987; 
Kuchuk et al., 1991) or the equivalent pressure point (Ozkan et ah, 1989; Daviau et al., 1988). 
In this work, the former method was used, with some adjustments due to the considered 
multilateral wellbore model. 

Each well ramification was segmented at 100 equidistant points. This number of well¬ 
bore segments was arbitrarily chosen. Then, for a given time, 100 pressure values were 
computed setting x = z = r w . In each layer, the term was evaluated applying, in y 
direction, the wellbore length corresponding to those 100 segments. Wellbore pressure was 
determined by taking the average between those 100 computed pressures. Figure 2 repre¬ 
sents an example where each wellbore ramification was divided into 4 different segments. 
The generalized expression for the average pressure in a wellbore divided into N s segments 
is given by: 


n b 
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Figure 2. 


Schematics showing the points used to compute the wellbore average pressure 


The constant flow-rates hypothesis, in its turn, is expected to be quickly satisfied (Gao, 
1987) and, therefore, should not yield significant errors. In fact, the application of equation 
(9) only requires the actual values of layer flow-rates to determine the overall mechanical 
skin factor. A simple but yet adequate estimate for the reservoir skin factor may be obtained 
replacing layer flow-rates by their respective flow capacities ( hj \Jk Xj k Vj ) in equation (9). 
Besides that, the only assumption made regarding layer flow-rates is that shall remain 
constant. 

3.4 Pressure derivative behavior 

Well testing in horizontal wells may present up to four distinct flow regimes (Goode & 
Thambynayagam, 1987; Rosa & Carvalho, 1989; Odeh & Babu, 1990; Kuchuk et al., 1991): 
early radial, intermediate linear, late radial and boundary governed flow. 

Early radial flow occurs at initial times, while the reservoir vertical boundaries are not 
felt by the wellbore. During this period, pressure behavior is equivalent to a fully penetrating 
vertical well in an infinite reservoir whose thickness is equal to the horizontal wellbore length. 
Pressure derivative with respect to the natural logarithm of time (hereafter referred to just 
as pressure derivative) is constant (Daviau et al., 1988). Based on the formulation for single¬ 
layer reservoirs (Ozkan et al., 1989), early time pressure derivative in a multilayer system is 
expected to stabilize at: 


1 q » ( 12 ) 

9In(f) 2 L t {k zx ) eq 

Since the considered reservoir model accounts for distinct layer properties, the pressure 
transient pulse propagates differently in each layer. Therefore, the intermediate linear flow 
might not simultaneously develop in all layers. Hence, the transition period between early 
and late radial flow is more accurately depicted by a 3-dimensional transitional flow rather 
than a linear flow regime per say. For this reason, pressure derivative presents no a priori 
characteristic signature. 

At late times, if the reservoir is sufficiently wide, a second radial flow regime develops 
(Kuchuk et al., 1991). Pressure response resembles a vertical well fully penetrating the 
reservoir thickness, inducing another constant pressure derivative level (Jelmert & Thomp¬ 
son, 1991; Kuchuk & Habashy, 1996). Pressure derivative, then, reflects the permeabilities 
in directions x and y : 


dAP _ 1 qp 
d ln(t) 2 ht{k' X y)eq 


(13) 


In equation (13), ht is the reservoir overall thickness and the ( k xy ) eq term denotes the 
thickness averaged permeability in the horizontal plane, defined as: 
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The last flow regime starts when the reservoir outer lateral boundaries are felt by the 
wellbore. In this work, a laterally infinite reservoir is considered. Thus, the boundary 
governed flow regime is never achieved, and the last identifiable flow regime is the late 
radial. 


4 Two-Phase Flow Through Multilateral Horizontal Wells 

The formulation for two-phase flow through horizontal wells in single-layer systems was 
reached by Peres and Reynolds (2003) and Peres et al. (2004). In this section, their solution 
will be introduced and, then, extended to multilayer stratified reservoirs. The main insight 
is to split the wellbottom hole pressure in two terms: one that encompasses the single-phase 
oil flow and another that derives from the mobility contrast between injected and underlying 
fluids (Barreto Jr. et al., 2011; Bela et al., 2019b). The reservoir is considered to be under 
the same assumptions made in section 3. Moreover, capillary effects are neglected. 

Similar to section 3, this section first introduces the existing solution for single-layer 
reservoirs with horizontal wells and then extends it to the multilayer case. 

4.1 Two-phase flow in single-layer reservoirs 

As proposed by Peres and Reynolds (2003), pressure during injectivity tests at a single¬ 
layer reservoir is computed as: 


A P(t) = A P und (t) + 


Qinj 


irh 


^und 
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min (dz,VF Z x) 


und | dv 

A t(r,t) ) r 


Xund ^\ dv 

Xt{r,t ) J r 


(15) 


In equation (15), A P un d{t) stands for the pressure change due to the single-phase 
displacement of the reservoir underlying fluid, determined as depicted in section 3.1. The 
symbols \ un d and A ( denote, respectively, the endpoint mobility of the reservoir underlying 
fluid and total mobility, as defined by Peres and Reynolds (2003) and Bela et al. (2019b). 
The terms between parenthesis represent the mobility contrast in the x direction, the xy 
plane, inside the damaged region and in the zx plane, respectively. The flood fronts in 
each plane are determined as depicted by Buckley and Leverett (1941), using the proper 
equations for each geometry. 

Equation (15) is reached under the assumption that flow-rate is constant within the 
flooded region. Peres and Reynolds (2003) stated the criteria required for this hypothesis 
to hold and showed that, in most practical conditions, the flood front is indeed within the 
steady-state region. Integration limits in equation (15) come from the flooding patterns 
proposed by Deppe (1961). The first integral upper limit is defined as: 
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{ 7T L 

(16) 

max(y,a:j?) 

For non centered wells, the effective flooded thickness varies in the x direction. That 
is the reason why thickness is expressed as a function of x in the first integral of equation 
(15). Figure 3 displays the flood front propagation model that was considered in this work. 
For values of x such that ^ < x < x aU x, the effective flooded zone thickness is estimated 
by a linear interpolation between the reservoir thickness and the smallest distance between 
the well and a vertical boundary (Peres et ah, 2004). The coordinate x aux is defined as: 


%aux — 




Tiy ln ( 


2 dz ) 


In 


2ndz sin(7r dz/h) 


l n ( _h_\ _ 1 


(17) 


It is interesting to notice that, as the flood front develops, the integrals in equation 
(15) become less or more significant, since the integral limits change in time (Deppe, 1961; 
Peres & Reynolds, 2003). 


4.2 Two-phase fluid injection in multilayer reservoirs 

Similar to what was done in section 3.2, the multilayer formulation begins by applying 
equation (15) at a given layer j: 


APj(t) — A P un dj[t) + 




k j A und-j h j 




dx 

hj{x) + 



Lj/2 




Considering a stratified reservoir plays a critical role at the proposed two-phase formu¬ 
lation. If there were formation crossflow, once the flood front at a given layer had reached 
a vertical boundary, it would possibly invade the adjacent layers. Flood front propagation, 
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0 g *ajx 



Figure 3. Linear flood front propagation from a non centered horizontal well 
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then, would be much more complicate to foresee. Since formation crossflow was not ac¬ 
counted for by the considered reservoir model, flood front in each layer may be computed 
independently and analogously to the single-layer model. 


Inspired by the analytical solution for vertical wells achieved by Barreto Jr. et al. (2011), 
a weighting variable Aj (t) is now introduced: 
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dx 
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(19) 


Replacing equation (19) in (18): 


A Pj(t) = A P undj (t) + q^AjA) 


( 20 ) 


As demonstrated in section 3.2, the term associated with the single-phase fluid dis¬ 
placement is, in fact, the same for all layers. For this reason, it will be hereafter denoted 
as A P und {t)- Besides, the considered model states that pressure change is also equal in all 
layers, except for the hydrostatic column. Thus, layer flow-rate may be computed as: 


= 


A P(t) - A P und (t) 
Aj it) 


Summing up all layer flow-rates: 
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Rearranging equation (22): 
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(23) 


Equation (23) depicts the pressure behavior in multilayer systems under fluid injection. 
In single-layer reservoirs, it yields the same result as equation (15), which assures that the 
suggested formulation is applicable to reservoirs with any number of layers. Furthermore, 
equations (21) and (23) are exactly the same that model two-phase flow in multilayer reser¬ 
voirs with vertical wellbores (Barreto Jr. et al., 2011). The only differences consist of the 
definition of Aj(t) and the computation of A P und (t). 

Equation (21) may be used to estimate layer flow-rate profile during the test. A com¬ 
putational implementation of the proposed analytical might use equation (21) to update 
layer flow-rates at each time step, using the following procedure: 


11 














manuscript submitted to 


1. Provide an estimate for layer flow-rates at the initial time step, for instance, using 
layer flow capacities (hj \Jk Xj k Vj ) 

2. Compute the single-phase underlying fluid contribution A P un d(t), the Aj(t) coeffi¬ 
cients and the wellbore pressure A P(t) using equations (9), (19) and (23) 

3. Update the time step 

4. Use equation (21) and the computed values of A P un d, Aj and A P from the previous 
time step to update flow-rates in all layers 

5. Repeat steps 2 to 4 

However, recalling the proposed formulation for single-phase flow developed in section 
3.2, equation (9) was reached assuming that layer flow-rates remain constant in time. This 
hypothesis, although not explicitly mentioned so far in this section, is also required to 
achieve equation (21), since the A P un ,d term is computed through equation (9). Therefore, 
the procedure described above to determine layer flow-rate profile is inconsistent with the 
model depicted in section 3.2, which assumes constant layer flow-rates. 

Nonetheless, whenever layer properties are remarkably different (Gao, 1987; Ehlig- 
Economides & Joseph, 1987), a rate-transient might develop. Even in reservoirs where layer 
properties are similar, flow-rates might change in time if the difference between layer skin 
factors is significant (Bela et al., 2019a). In section 5.3, an analysis is conducted over the 
influence of assuming constant layer flow-rates. 

4.3 Pressure derivative behavior 

Peres and Reynolds (2003) showed that two-phase flow in reservoirs with horizontal 
wells might present several distinct flow regimes, which depend on the development of the 
flooded region and the pressure transient zone. 

They denoted the first regime as first radial-first radial, since both the pressure transient 
pulse and the flood front are propagating in the vertical plane perpendicular to the wellbore. 
During this flow regime, pressure derivative reflects the underlying fluid properties, and 
behaves similar to the early radial regime observed at a single-phase flow. Therefore, pressure 
derivative during this flow regime attains a constant level that is related to the equivalent 
properties in the vertical plane: 


dAP _ 1 q 

Shi(i) ~ 2 L t (k zx ) eq X U nd 


(24) 


After the early radial flow, a linear regime develops in each layer. However, as mentioned 
in section 3.4, the linear flow might start and end at different times in each layer. Therefore, 
it is not possible to state a priori that the characteristic derivative signature related to this 
flow regime will be identified. 

Again, the linear regime is followed by a late radial flow. Since the reservoir is assumed 
to be laterally infinite, eventually pressure diffusion will simultaneously propagate through 
the horizontal plane in all layers. Pressure derivative during this late radial flow period will 
depend on the flood front radius in each layer. Using the same notation from Peres and 
Reynolds (2003), either a ’’second radial - first radial” or a ’’second radial - first linear” or 
a ’’second radial - second radial” flow regime may be observed. Assuming that the ’’second 
radial - first radial” flow develops in all layers, pressure derivative is approximated by: 


dAP ~ 1 q ( | ht (1 -M) \ 
dhi(i) ~ 2 h t (k xy ) eq \ und \ U. M ) 
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where M stands for the endpoint mobility ratio. The expression for pressure derivative 
while the ’’second radial - second radial” regime occurs in all layers is given by: 


dAP _ 1 q / | ht (1 -M) \ 
d\n(t) ~ 2 ht ( kxv ) eq X und \ 2 L t M M J 


Approximations for the other flow regimes may be obtained by adapting the generalized 
expression for pressure derivative in single-layer reservoirs achieved by Peres and Reynolds 
(2003). 

5 Results and Discussion 

To assess the accuracy of the formulation detailed in sections 3 and 4, a set of cases 
was run on a finite-differences-based numerical simulator and results were compared to the 
analytical data. 

The grid considers the reservoir width (x-direction) equals 5.0 km (3.11 mi), whereas 
the reservoir length (y-direction) equals 2.5 km (1.55 mi). To mitigate spatial discretization 
effects a 5x5x5 cartesian refinement was applied around the wellbore. It was assumed 
that the well is open during a 96h period. In all cases, a 0.108 m wellbore radius was 
employed. Values for the underlying fluid and rock compressibilities were 1.14-10 -4 cm 2 /kgf 
and 8.0-10 -5 cm 2 /kgf, respectively. 

In all cases, it was considered that well drilling and completion change the permeability 
around it, representing the skin effect. Grid blocks containing the wellbore were built with 
0.6 m width and 0.6 m height (that is, distance between the wellbore and a grid block face 
is 0.3m). Permeability inside those blocks was changed to represent the damaged region. 
Skin zone in the analytical solution was modeled defining the skin radius required so that 
the cross-sectional damaged region would present the same area as the numerical simulator, 
yielding a skin radius of 0.35 m. Figure 4 illustrates this procedure. 

Table 1 shows the reservoir parameters employed at the tested cases. Cases Al to A4 
represent the same reservoir, but cases A2 and A4 account for vertical anisotropy. Cases 
Al, A3, B and C, in their turn, consist of isotropic systems. 

5.1 Single-phase flow results 

The results associated with cases Al, A2, B and C are displayed at figures 5 and 6. 
The single-phase results and analysis for cases A3 and A4 were similar to cases Al and A2, 
so they were omitted to avoid redundancy. Analytical model data is represented by solid 
lines, while the dots stand for the numerical data. Theoretical constant pressure derivative 



Damaged region 


Figure 4. Damaged Zone Representation 
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Case 

q 


<t>j 

k X yj 

k Z j 

Sj 

dzj 

Lj 

h.j 

(m 3 /d) 

(cP) 

(mD) 

(mD) 

(m) 

(m) 

(m) 

A1 

1000 

4.5 

0.18 

500 

500 

3.5 

7.5 

50 

15 

0.18 

750 

750 

3.5 

10.0 

150 

20 

A2 

1000 

4.5 

0.18 

500 

130 

4.7 

7.5 

50 

15 

0.18 

750 

190 

7.6 

10.0 

150 

20 

A3 

1000 

1.3 

0.18 

500 

500 

3.5 

7.5 

50 

15 

0.18 

750 

750 

3.5 

10.0 

150 

20 

A4 

1000 

1.3 

0.18 

500 

130 

4.7 

7.5 

50 

15 

0.18 

750 

190 

7.6 

10.0 

150 

20 

B 

1500 

5.7 

0.23 

500 

500 

3.5 

10.0 

100 

20 

0.27 

1450 

1450 

3.5 

5.0 

120 

15 




0.28 

450 

450 

3.5 

7.5 

120 

15 

C 

1200 

0.9 

0.15 

680 

680 

3.5 

10.0 

120 

20 




0.22 

830 

830 

3.5 

5.0 

120 

10 


Table 1 . Reservoir parameters 


levels associated with the early and late radial flow regimes were determined from equations 
(12) and (13). These levels are plotted in dark dashed lines. 

Despite some divergences, mainly at early times, all cases presented a close agreement 
between the analytical model and the numerical simulator. In cases A1 and A2, pressure 
derivative clearly attains two distinct constant levels, matching the corresponding theoretical 
values foreseen by equations (12) and (13). In cases B and C, the early radial regime is 
noticeably shorter. This is an immediate consequence of layer properties, since the most 
permeable layer in cases B and C presents shorter thickness and wellbore length, compared 
to the most permeable layer in cases A1 and A2. Therefore, in cases B and C, the reservoir 
vertical boundaries are felt earlier. In fact, pressure derivative profile indicates that, in case 
C, early radial flow lasts only a few seconds. 




Figure 5. Single-phase pressure and pressure derivative profile for cases A1 (left) and A2 (right) 
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Figure 6. Single-phase pressure and pressure derivative profile for cases B (left) and C (right) 


The transitional regime between early and late radial flows resembles a linear flow, as 
observed in single-layer reservoirs with horizontal wells. However, as mentioned in section 
3.3, it is not possible to assure that this flow regime actually develops in multilayer reservoirs. 
The main reason for that is the difference between layer properties. Thereby, the linear flow 
regime might not simultaneously develop in all layers. 

As depicted in section 3.2, the proposed formulation considers that layer flow-rate do 
not change during the test. Therefore, the analytical model implementation assumed that 
layer flow-rates would remain constant. The numerical simulator, in its turn, allows the 
existence of a rate transient. The formulation for vertical wellbores suggests that layer flow- 
rates should stabilize after a short transient period (Gao, 1987; Elrlig-Economides & Joseph, 
1987). Thus, the analytical model is expected to present more relevant errors at early times 
and a higher accuracy during the late radial flow. 

Figures 7 and 8 portray the flow-rate profiles computed by the numerical simulator, as 
well as the theoretical flow-rate levels according to each layer flow capacity. In all cases, layer 
flow-rates change at early times. This behavior is expected, due to the differences between 
layer properties. Moreover, the rate transient is likely related to the distinct identifiable 
flow regimes (early radial, 3-dimensional transitional regime and late radial). 




Figure 7. Single-phase flow-rate profiles of cases A1 (left) and A2 (right) 
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t(h) 


t(h) 


Figure 8. Single-phase flow-rate profiles of cases B (left) and C (right) 


Nonetheless, in all cases, layer flow-rates stabilized after this short transient period. 
Although the proposed single-phase formulation only requires that layer flow-rates remain 
constant, no assumption is made about which level the flow-rates should attain. In fact, 
the rate profile displayed in figures 7 and 8 might suggest that the analytical model would 
present significant divergences at early times, while layer flow-rates are changing. Regard¬ 
less, numerical and analytical data were close even during the rate transient period. 

The main divergences between numerical and analytical data observed at early times 
are probably related to the simulation grid instead. Numerical simulation is sensitive to the 
employed grid refinement. Besides, the damaged region in the simulation grid is modeled by 
cubes, although the skin zone is in fact cylindrical. These issues are more likely the reason 
for the observed early time mismatch, rather than an eventual inaccuracy of the analytical 
model. 

5.2 Equivalent permeability estimates from single-phase flow data 

As described in section 3.4, two equivalent permeabilities might be estimated from the 
two constant pressure derivative levels. Those equivalent permeabilities are not necessarily 
equal, even for the isotropic cases, since they reflect differently weighted averages, as shown 
in equations (12) and (13). Thereby, the constant pressure derivative levels reached by the 
suggested analytical solution were applied to compute the reservoir equivalent permeability. 

The estimated equivalent permeabilities for cases Al to C might be found in table 2. In 
all cases, the reservoir equivalent permeability could be determined with a good accuracy, 


Case 

( k ; 

Real 

zx)eq 

Est. 

(mD) 

Error (%) 

( kxy)eq 
Real Est. 

(mD) 

Error (%) 

Al 

688 

717 

4.2 

643 

631 

-1.8 

A2 

347 

366 

5.6 

643 

628 

-2.4 

B 

1018 

1057 

3.8 

907 

912 

0.5 

C 

653 

— 

— 

637 

632 

-0.7 


Table 2. Computed equivalent permeabilities - single-phase cases 
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in both radial flow regimes. In case C, the equivalent permeability associated with the early 
radial flow could not be estimated, due to the extremely short duration of this flow regime 
in case C. 

As previously commented, the proposed formulation assumes constant layer flow-rates, 
although a rate transient might occur at early times. For this reason, estimated values 
for the equivalent permeability in the z ~ x plane were expected to present higher errors. 
Nonetheless, ( k zx ) eq could be computed with good precision in all cases where the early 
radial flow lasted long enough to yield a clear constant pressure derivative portion. On the 
other hand, equivalent permeabilities related to the late radial regime were determined with 
very good accuracy in all cases. 

In a real field test, however, the constant pressure derivative level related with early ra¬ 
dial flow might not be identifiable in reservoirs where the wellbore storage effect is relevant. 
Thereby, the proposed formulation could provide close estimates for the equivalent perme¬ 
abilities associated with the flow regime that is most important for practical applications. 

5.3 Two-phase flow results 

The developed solution for two-phase flow was assessed considering that the injected 
fluid viscosity equals 0.52 cP in all cases. It was also assumed that all layers present the 
same relative permeability curves. These curves are exhibited in figure 9. Applying the 
criteria suggested by Peres and Reynolds (2003), it is possible to check that the injection 
flow-rate is low enough to assure that the flood front is always within the steady state region 
in all layers of all cases. 

Two computational implementations of the proposed two-phase model were developed. 
The first considered that layer flow-rates would remain constant during the test, whereas 
the second updates layer flow-rates at each time step, according to the procedure detailed 
in section 4.2. Figure 10 portrays the results for pressure, pressure derivative and layer 
flow-rates using both implementations. Reservoir properties are the same of case Al. Two 
constant flow-rate levels were also plotted. Those levels were determined according to layer 
flow capacities {h 3 ~Jk Xj k y .) and from the product between permeability in zx plane and 
wellbore length (L ? ^k Zj k Xj ). Based on the pressure derivative profile depicted in section 3.4, 
such flow-rate levels are likely attained during late radial and early radial flow, respectively. 

Figure 10 evidences that differences between both implementations is almost unde¬ 
tectable. Similarity between the two implementations is consistent with the analytical model 



Figure 9. Relative permeability data 
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Figure 10. Comparison between the different two-phase model implementations considering 
reservoir properties from case A1 


developed in sections 3.2 and 4.2. The actual values of layer flow-rates are only required 
to determine the reservoir overall skin factor in equation (9). Other than that, the pro¬ 
posed solution do not demand explicit knowledge of individual layer flow-rates. Moreover, 
it is interesting to highlight that, although the transient rate implementation allows layer 
flow-rates to change in time, they remained approximately constant during the entire test. 

Similar results (not shown here) were obtained for the other cases. Thereby, only the 
constant flow-rate implementation was compared with the numerical simulator. The reason 
for choosing this implementation over the transient rate one is because it is more consistent 
with the suggested single-phase flow model. 

The two-phase flow results for cases Al to C are displayed in figures 11 to 13, as well 
as the theoretical pressure derivative levels. Early radial derivative level was determined 
from equation (24). During late radial flow, it was assumed that the ’’second radial - first 
radial” flow regime occurs in all layers. Thus, the constant derivative level was computed 
from equation (25). 




Figure 11. Two-phase pressure and pressure derivative profile for cases Al (left) and A2 (right) 
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Figure 12. Two-phase pressure and pressure derivative profile for cases A3 (left) and A4 (right) 


A constant pressure derivative associated with late radial flow is indeed observed in all 
cases. However, flood front profiles (not displayed here) show that cases Al and A2 were the 
only cases where the ’’second radial - first radial” regime lasts until the end of the test. In 
these cases, the reservoir vertical boundaries are not reached by the flood front by the end 
of the test. In all other cases, a ’’second radial - first linear” flow starts at the narrower layer 
around 50 to 60 h, whereas the thicker layer remains at the ’’second radial - first radial” 
regime. For this reason, pressure derivative curve detaches from the theoretical pressure 
derivative level at its final portion, as a new transitional flow regime develops. 

Early radial flow, on the other hand, is not clearly identified in any case. Moreover, 
significant mismatches between analytical and numerical data were, once again, observed 
at early times. Such divergences are explained by the same reasons detailed in section 5.1. 
Furthermore, pressure response for the two-phase flow is even more sensitive to the employed 
simulation grid, since the flood front computation gets more accurate as the simulation grid 
is more refined. The mobility contrast between injected and underlying fluids also helps to 
understand why the early constant derivative level is not attained by pressure derivative, 
since the injected fluid is less viscous than the reservoir underlying fluid. 




Figure 13. Two-phase pressure and pressure derivative profile for cases B (left) and C (right) 
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Figure 14. Two-phase flow-rate profiles of cases A1 (left) and A2 (right) 


In reservoirs subject to two-phase flow, a characteristic pressure derivative signature 
is noticed, as the flood front overcomes the damaged region. Such signature consists of a 
sharp shift in pressure derivative, which might attain negative values (Peres & Reynolds, 
2003; Bela et al., 2019a). This feature is easily noticed in all cases. However, the most 
pronounced divergences between numerical and analytical models occurred precisely while 
the flood front was crossing the skin zone radius. In fact, pressure curves clearly detach 
from one another while the blunt derivative shift is observed. 

Once again, the main reason for that relies on the skin zone representation at the 
simulation grid representation. Whereas the analytical solution accounts for a cylindrical 
damaged region, it is modeled by cubic grid cells in the numerical simulator. Pressure 
derivative profile highlights the differences between damaged zone representations. Injec¬ 
tivity tests in multilayer reservoirs with formation damage might present several derivative 
shifts, as the flood front overcomes the skin radius in each layer (Bela et al., 2019a). Al¬ 
though the numerical derivatives from cases A3, A4 and C present only one sharp drop, 
several shifts might be identified at the analytical derivative curves. The proposed solution 
also presents an evident derivative shifts in cases Al and B, while numerical data, quite 
unexpectedly, shows no derivative bump. 




Figure 15. Two-phase flow-rate profiles of cases A3 (left) and A4 (right) 
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Figure 16. Two-phase flow-rate profiles of cases B (left) and C (right) 


As previously mentioned, the suggested formulation is achieved by assuming constant 
layer flow-rates, but the numerical simulator accounts for rate transient. Simulated layer 
flow-rate profiles, portrayed in figures 14 to 16, evidences that the rate transient period 
is indeed quite short, similar to the single-phase results (figs 7 and 8). The theoretical 
flow-rate levels were again determined according to layer flow capacities. Early time rate 
transient may also explain the analytical and simulated data mismatch at initial times, in 
addition to the other reasons mentioned above. However, in all cases, layer flow-rates quickly 
stabilized. Furthermore, comparison between figures 10 and 14 shows that layer flow-rate 
history was not replicated by the transient rate analytical implementation, that uses the 
procedure depicted in section 4.2. 

5.4 Equivalent permeability estimates from two-phase flow data 

Constant pressure derivative levels obtained with the proposed analytical model might 
be used to estimate the reservoir equivalent permeabilities from equations (24) to (26). Since 
the constant derivative level related to early radial flow was not clearly identified in any case, 
no estimate was computed for ( k zx ) eq . As commented in section 5.2, early radial flow may 
be concealed by wellbore storage effect in a real field test. Furthermore, this flow regime 
is more sensitive to measurement noise. Thereby, a practical application of the suggested 
formulation may not be able to provide an estimate for ( k zx ) eq due to this reasons as well. 


Case 


( J^xyjeq 

Error (%) 

Real 

Est. 

A1 

643 

634 

-1.5 

A2 

643 

678 

5.4 

A3 

643 

632 

-1.8 

A4 

643 

609 

-5.2 

B 

907 

1061 

16.9 

C 

637 

643 

1.0 


Table 3. Computed equivalent permeabilities - two-phase cases 
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Equivalent permeability in xy plane, in its turn, was estimated using equation (25), 
since the late time flow regime was classified as ’’second radial - first radial” according to 
the flood front profiles. Thus, in cases A3, A4, B and C, only the constant derivative portion 
corresponding to the ’’second radial - first radial” was employed to estimate the equivalent 
permeability (that is, the final derivative portion after the derivative detaches from the 
constant level was discarded). 

The computed equivalent permeabilities are exhibited in table 3. In case B, late radial 
flow was shorter, due to the beginning of a linear flood front propagation inside the narrower 
layer. Thus, pressure derivative did not attain the theoretical level forecast by equation 
(25), and equivalent permeability estimate was compromised. However, to overcome this 
problem a longer testing time may be employed. Then, pressure derivative related with 
the ’’second radial - second radial” flow regime may be employed to estimate the reservoir 
equivalent permeability. Results for the other cases presented a good accuracy, stating that 
the proposed solution may provide helpful information about the reservoir. 


6 Conclusions 

A real space solution for single-phase flow in multilayered reservoirs with multilateral 
horizontal wells was achieved. The proposed formulation is based on the real space solutions 
for single-layer reservoirs developed by Daviau et al. (1988); Ozkan et al. (1989). The 
suggested solution is achieved assuming a uniform influx wellbore. An adjustment was made 
at the average pressure technique, so that the developed formulation could be applied to 
compute pressure response from an infinite conductivity wellbore (which is a more realistic 
hypothesis than uniform influx). 

The proposed single-phase model was extended to depict two-phase flow inside a reser¬ 
voir subject to the same assumptions. Inspired by the single-layer formulation reached by 
Peres and Reynolds (2003), the proposed solution splits the measured pressure change in 
two terms: one related to the mobility differences between injected and underlying fluids, 
and another related to the single-phase flow that occurs beyond the flood front. To the best 
of the authors’ knowledge, the developed real space single-phase formulation has never been 
presented in literature before, whereas the achieved two-phase model is the first solution 
ever proposed for two-phase flow in multilayer reservoirs with multilateral wells. 

The suggested formulation accuracy was assessed by comparison with a commercial 
numerical flow simulator. Results for a set of cases showed a good agreement both for the 
single-phase and two-phase models. Some early time divergences were noticed, particularly 
while the flood front overcomes the damaged regions in the two-phase cases. Nevertheless, 
the overall pressure behavior showed a close agreement between analytical and numerical 
data. These results not only assure that the suggested formulation is valid, but also prove 
the efficiency of the proposed adjustments to the average pressure technique. Furthermore, 
the developed solution was applied to estimate the reservoir equivalent permeability with 
good accuracy in all cases where a constant pressure derivative level was noticed. 
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Nomenclature 

A Weighting variable 

b Upper integral limit in equations (15) and (18) 
c Compressibility 

dz Smallest distance between the wellbore and a vertical boundary 
erf Error function 

exp Exponential function 
h Thickness 

k Permeability 

L Wellbore length 

In Natural logarithm function 

M Endpoint mobility ratio 

P Pressure 

q Flow-rate 

r Radius 

S Skin factor 

sin Sine function 

t Time 

Hydraulic diffusivity coefficient 
A Fluid mobility 

p Viscosity 

</> Porosity 

a Auxiliary parameter 
r Silent time integration variable 


Subscripts 

eq Equivalent 
F Related to the flood front 
i Initial 

inj Related to the injected fluid 
j Related to layer j 
skin Related to the damaged region 
t Total 

und Related to the reservoir underlying fluid 
x With respect to x direction 
xy With respect to xy plane 
y With respect to y direction 
0 With respect to z direction 
zx With respect to zx plane 
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